clc,clear,close all;
%%
% t = 0:0.01:10;
% y = 3*sin(2*pi*10*t)+sin(2*pi*40*t)+sin(2*pi*200*t)+5;
% 
% %plot(t,y,'b-');
% 
% 
% axis([0 1 0 10]);
% 
% datalength = length(y);
% 
% yFFT = fft(y,datalength);
% 
% mag = abs(yFFT/datalength);
% 
% magy= mag(1:datalength/2+1);
% 
% figure;
% 
% plot(magy,)
% 生成示例信号（含50Hz和120Hz正弦波+噪声）
% Fs = 1000;            % 采样频率1000Hz
% t = 0:1/Fs:1;         % 1秒时间向量（1001个点）
% f1 = 50; f2 = 120;    % 信号频率成分
% %signal = 0.7*sin(2*pi*f1*t) + sin(2*pi*f2*t); 
% %noise = 2*randn(size(t)); % 高斯噪声
% data = 3*sin(2*pi*10*t)+sin(2*pi*40*t)+sin(2*pi*200*t)+5; % 含噪信号
% 
% % 傅里叶变换
% N = length(data);             % 信号长度
% Y = fft(data);                % 执行FFT
% P2 = abs(Y/N);                % 双侧频谱
% P1 = P2(1:N/2+1);             % 取单侧频谱（0~Nyquist频率）
% P1(2:end-1) = 2*P1(2:end-1);  % 修正幅度（除直流分量外）
% 
% % 创建频率轴
% f = Fs*(0:(N/2))/N; 
% 
% % 绘制频谱图
% figure;
% plot(f, P1, 'LineWidth', 1.5);
% title('单边幅度谱');
% xlabel('频率 (Hz)');
% ylabel('幅度');
% grid on;
% xlim([0 200]);  % 聚焦0-200Hz范围
%%
% t = 0:0.01:10;
% y = 3*sin(2*pi*10*t)+sin(2*pi*40*t)+sin(2*pi*200*t);
% %数据长度
% datalength = length(y);
% %对数据进行快速傅里叶变换
% yFFT = fft(y,datalength);
% mag = abs(yFFT);
% %幅值归一化
% mag = mag*2/datalength;
% pha = angle(yFFT)*180/pi;
% %无效相位 置0
% for i = 1:datalength
%     if (mag(1,i)<0.3)
%         pha(1,i) = 0;
%     end
% end
% %序列坐标
% n = 0:datalength-1;
% %采样频率
% fs = 1/0.01;
% %序列频率
% f = (0:datalength-1)*fs/datalength;
% 
% figure;
% stem(f(1:datalength/2),pha(1:datalength/2),'b-','LineWidth');


%%
%信号的合成与分解
%dot(a,y)是向量a和向量y的内积
%dot(a,y)/dot(y,y)成为向量a在向量y上的投影
% a = [80 90 70];
% x = [1 0 0 ];
% y = [0 cos(pi/3) sin(pi/3)];
% z = [0 -sin(pi/3) cos(pi/3)];
% x_coor = sum(a.*x)/sum(x.*x);
% y_coor = dot(a,y)/dot(y,y);
% z_coor = dot(a,z)/dot(z,z);

%% 
% 方波的合成与分解
% 函数的内积：两个函数所在区间的积分。
% 函数正交：内积为0，两个函数正交。
% 正交函数集：类比于正交基。
% 三角函数集的正交性
% {1，cosx，sinx，cos2x,sin2x,...,cosnx,sinnx}内的函数在区间[-pi,pi]上彼此正交。
% 周期信号傅里叶级数
% 周期信号f(t),其周期为T，但满足狄里赫利条件时，它可以分解为如下三角级数，该三角级数成为傅里叶级数。

%绘图变量
font_size = 12;axis_size = 10;line_width = 2;legend_size = 8;marker_size = 6;
figure_width = 14;figure_height = 8;BiaValue = 0;

%绘制方波信号
%基波频率
T0 = 4;w0 = (2*pi)/T0;
%信号值为1的区间
T1 = 1;
%绘图的周期个数（2*n+1）
n = 1.5; 

%绘图时长
t = (-n*T0-T0/2):0.001:(n*T0+T0/2);
%方波信号
y = zeros(1,length(t));
for i = 1:length(t);
    %将周期变化归一化
    t_ = double(abs(t(1,i)))-double(int32(abs(t(1,i)))/int32(T0))*double(T0);
    
    if(abs(t_)<T1)
        y(1,i) = 1;
    else
        y(1,i) = 0;
    end
end

%% 绘图
figure
set(gcf,'Unit','centimeters','Position',[3 3 3+figure_width 3+figure_height]);
set(gca,'LooseInset',get(gca,'TightInset')+[BiaValue,0,0,0],'FontName','Times New Roman','FontSize',font_size);
plot(t, y, 'b-', 'LineWidth', line_width)
set(gca,'FontSize',font_size,'FontName','Times New Roman')
xlim([(-n * T0 - T0 / 2) (n * T0 + T0 / 2)])
ylim([0 2])
set(gca,'XTick',(-n * T0 - T0 / 2) : 4 : (n * T0 + T0 / 2))
set(gca,'YTick',0 : 1 : 2)
xlabel('时间 (s)', 'FontSize', font_size, 'FontName', 'Times New Roman')
ylabel('信号幅值', 'FontSize', font_size, 'FontName', 'Times New Roman')
title('方波信号波形', 'FontSize', font_size*1.5, 'FontName', 'Times New Roman')
 
%% 方波信号的合成
% 直流分量系数
c0 = 2 * T1 / T0;
% 拟合的谐波个数
num = 7;
% 谐波分量对应的系数
ck = zeros(1, num);
for k = 1:num
    ck(1, k) = sin(k * w0 * T1)/ (k * pi);
end
 
% 绘制单边谱（将复指数分解转换为余弦信号分解）
ak_cos = 2 * ck;
 
%% 绘图：绘制单边谱
figure
set(gcf,'Unit','centimeters','Position',[3 3 3+figure_width 3+figure_height]);
set(gca,'LooseInset',get(gca,'TightInset')+[BiaValue,0,0,0],'FontName','Times New Roman','FontSize',font_size);
stem(0:1:num, [c0 ak_cos], 'b-', 'LineWidth', line_width)
set(gca,'FontSize',font_size,'FontName','Times New Roman')
xlim([0 num])
ylim([-0.25 0.75])
set(gca,'XTick',0 : floor(num/3) : num)
set(gca,'YTick',-0.25:0.25:0.75)
xlabel('频率 (单位间隔为基波频率)', 'FontSize', font_size, 'FontName', 'Times New Roman')
ylabel('幅值', 'FontSize', font_size, 'FontName', 'Times New Roman')
title('方波信号的单边幅度谱', 'FontSize', font_size*1.5, 'FontName', 'Times New Roman')
 
%% 绘图：演示分解出的信号
figure
set(gcf,'Unit','centimeters','Position',[3 3 3+figure_width*1.2 3+figure_height]);
set(gca,'LooseInset',get(gca,'TightInset')+[BiaValue,0,0,0],'FontName','Times New Roman','FontSize',font_size);
% 绘制直流分量
plot(t, c0*ones(1,length(t)), 'k-', 'LineWidth', line_width)
hold on
% 绘制谐波分量
% 自动调整颜色
greylevel = 0.1:(0.9-0.1)/(num-1):0.9;
for k = 1:num
    plot(t, ak_cos(1, k)*cos(k*w0*t), 'Color', [greylevel(1,k) greylevel(1,k) greylevel(1,k)], 'LineStyle', '-.', 'LineWidth', line_width)
end
set(gca,'FontSize',font_size,'FontName','Times New Roman')
xlim([(-n * T0 - T0 / 2) (n * T0 + T0 / 2)])
ylim([-1 1])
set(gca,'XTick',(-n * T0 - T0 / 2) : 4 : (n * T0 + T0 / 2))
set(gca,'YTick',-1:0.5:1)
xlabel('时间 (s)', 'FontSize', font_size, 'FontName', 'Times New Roman')
ylabel('信号幅值', 'FontSize', font_size, 'FontName', 'Times New Roman')
% 增加图例
str = [];
str{1} = sprintf('直流分量');
for i = 1:num
    str{i+1} = sprintf('%d次谐波分量', i);
end
legend(str, 'FontSize', legend_size, 'FontName', 'Times New Roman', 'Location', 'NorthEastOutside');
title('方波信号的直流及谐波分量', 'FontSize', font_size*1.5, 'FontName', 'Times New Roman')
 
%% 绘图：演示合成效果
figure
set(gcf,'Unit','centimeters','Position',[3 3 3+figure_width 3+figure_height]);
set(gca,'LooseInset',get(gca,'TightInset')+[BiaValue,0,0,0],'FontName','Times New Roman','FontSize',font_size);
% 原信号
plot(t, y, 'b-', 'LineWidth', line_width)
hold on
% 合成的信号
y_sum = c0*ones(1,length(t));
for i = 1:length(t)
    for k = 1:num
        y_sum(1,i) = y_sum(1, i) + ak_cos(1, k)*cos(k*w0*t(1, i));
    end
end
plot(t, y_sum, 'r-.', 'LineWidth', line_width)
set(gca,'FontSize',font_size,'FontName','Times New Roman')
xlim([(-n * T0 - T0 / 2) (n * T0 + T0 / 2)])
ylim([-1 2])
set(gca,'XTick',(-n * T0 - T0 / 2) : 4 : (n * T0 + T0 / 2))
set(gca,'YTick',-1:1:2)
xlabel('时间 (s)', 'FontSize', font_size, 'FontName', 'Times New Roman')
ylabel('信号幅值', 'FontSize', font_size, 'FontName', 'Times New Roman')
title('方波信号及合成信号', 'FontSize', font_size*1.5, 'FontName', 'Times New Roman')


